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Abstract. We present a simple, accurate method for computing singular or nearly sin¬ 
gular integrals on a smooth, closed surface, such as layer potentials for harmonic func¬ 
tions evaluated at points on or near the surface. The integral is computed with a reg¬ 
ularized kernel and corrections are added for regularization and discretization, which 
are found from analysis near the singular point. The surface integrals are computed 
from a new quadrature rule using surface points which project onto grid points in coor¬ 
dinate planes. The method does not require coordinate charts on the surface or special 
treatment of the singularity other than the corrections. The accuracy is about 0(/t 3 ), 
where h is the spacing in the background grid, uniformly with respect to the point of 
evaluation, on or near the surface. Improved accuracy is obtained for points on the 
surface. The treecode of Duan and Krasny for Ewald summation is used to perform 
sums. Numerical examples are presented with a variety of surfaces. 

AMS subject classifications: 65R20,65D30,31B10,35J08 

Key words: boundary integral method, layer potential, nearly singular integrals, Laplace equa¬ 
tion, surface integral, implicit surface 


1 Introduction 

We present a simple, accurate method for computing singular or nearly singular integrals 
defined on a smooth, closed surface in three-space. This method can be used to evaluate 
single or double layer potentials for harmonic functions, or the velocity and pressure in 
Stokes flow due to forces on a surface. The point of evaluation could be on or near the 
surface. To evaluate the integral, the kernel is first replaced by a regularized version. 
A preliminary value is found using a new quadrature rule for surface integrals which 
has the advantage that it does not require coordinate systems or a triangulation on the 
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surface. Instead we sum values of the integrand over quadrature points which project 
onto grid points in coordinate planes, in a way that would be high order accurate if the 
integrand were smooth. Corrections for the regularization and discretization are then 
added to achieve higher accuracy. These corrections are given by explicit formulas de¬ 
rived using asymptotic analysis near the singularity as in Q. The resulting value of the 
integral has 0(h p ) accuracy for p < 3, uniformly for points of evaluation near the surface, 
where h is the grid spacing in R 3 . For points on the surface, the accuracy is significantly 
improved by using a special regularization; see Section 3.3. For efficient summation we 
use the treecode algorithm of Duan and Krasny [111 designed for kernels with Gaussian 
regularization. The method presented here could be used, for example, to find values of 
the potential at grid points in R 3 close to the surface. It should be applicable to compu¬ 
tations with moving surfaces for which good accuracy is needed without extensive work 
to represent the surface at each time step. Other kernels could be treated by the same ap¬ 
proach, and the more accurate version of the method for computing values on the surface 
could be used for a variety of problems which can be formulated as integral equations. 

The present approach is an improvement and extension of the grid-based boundary 
integral method of Q. In the earlier work the integral was replaced by sums in coordi¬ 
nate charts using a partition of unity. The need for explicit coordinate systems requires 
knowledge of the surface that might be difficult to obtain for a moving surface. Further¬ 
more, if the coordinate system is too distorted, the accuracy will be poor because the 
discretization error will fail to be controlled by the regularization. Here we avoid these 
disadvantages by using a more direct rule for computing surface integrals, which was 
introduced for smooth integrands in [30[. This quadrature rule uses projections on coor¬ 
dinate planes rather than coordinate charts. Given a rectangular grid in R 3 , quadrature 
points are chosen as those points on the surface which project onto grid points in the co¬ 
ordinate planes, for which the normal to the surface has direction away from the plane. 
Weights for the quadrature points are found from a partition of unity on the unit sphere, 
applied to the normal vector at the point. The weight functions on the sphere are cho¬ 
sen universally, and do not depend on the particular surface. The resulting quadrature 
rule for surface integrals has high order accuracy, as allowed by the smoothness of the 
integrand and the surface. In effect, the method uses the existence of coordinate patches 
without having to refer to them explicitly. The quadrature points can be found efficiently 
if, for example, the surface is given analytically or numerically as the level set of a func¬ 
tion. Examples in |30| with smooth integrands illustrate the accuracy of this method with 
a variety of surfaces, including ones of large genus. 

For a variety of problems in partial differential equations, solutions can be written as 
integrals over surfaces, using a known fundamental solution. These include harmonic 
functions, electromagnetic waves, and viscous fluid flow modeled by the Stokes equa¬ 
tions. The specific representation makes this formulation an attractive approach for nu¬ 
merical methods. Often an integral equation on the surface must be solved, and much 
attention has focused on such problems. For the most familiar case, an integral equa¬ 
tion for the Dirichlet problem for the Laplacian, using the double layer potential, it was 
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proved in |4| that the discrete version of the integral equation using the present method 
has a unique solution and that it converges to the exact solution under grid refinement. 
The evaluation of the integral at points near the surface, in contrast to evaluation on the 
surface, generally requires extra care and can also be important in applications. If val¬ 
ues of a potential are needed at grid points in the computational domain, we can first 
compute values at points near the surface directly as nearly singular integrals. Once this 
is done, values at other grid points can be found in a cheap way by inverting a discrete 
Laplacian, as in pO] , In this paper the values obtained at the irregular grid points near 
the surface using integration have accuracy about 0(/j 3 ), while those found at the regular 
points are about 0(h 2 ) accurate; more accurate values could be found with more work. 


The most widely used numerical technique for integral formulations is the boundary 
element method (e.g. |l 14 24 25]|); the boundary is triangulated and matrix elements are 
computed for the integral operator on the surface using special quadrature rules. This 
method is especially useful for electromagnetic problems in which the surface does not 
evolve and may have corners and edges. Direct quadrature based on triangulations is 
also used, and can be accurate in the (exactly) singular case for smooth surfaces. Such 
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methods have long been used for modeling in chemistry and biology [ 23.24,361 
methods use corrected quadrature weights p9| or analytical evaluation of the singular 
part, treating the remainder in a standard way [16.30]. A careful direct quadrature (or 
Nystrom) method, introduced in [7j for electromagnetics, and further developed in [32], 
uses a partition of unity to reduce the integral to coordinate patches; a special patch in 
polar coordinates is used near the singularity. A different approach [131 is a spectral 
method, using spherical harmonics, assuming the surface can be mapped to a sphere. 
This method, applied in [ 12 29], appears to be advantageous when there are many bound¬ 
aries that are not greatly deformed. For use with many vesicles, the authors of [29 ] report 
difficulties in keeping the method of [7.32 [ accurate, perhaps because of the cut-off func¬ 
tion needed for the special patch near the singularity. Careful integral methods have the 
promise of calculating integrals very accurately even for complicated surfaces lacking 
smoothness; see |27) . An alternative approach, the kernel-free boundary integral 

method [34,35], replaces the calculation of the integral by the solution of an interface 
problem on a regular grid. This approach has the advantage that explicit knowledge 
of the integral kernel is not needed, and thus it can be applied to more general partial 
differential equations. 

One advantage of the present method is its simplicity. Detailed information about 
the surface is not required. No special treatment is needed near the singularity, except 
for the corrections which are added after the summation, using analytical formulas. By 
design, the errors are uniform with respect to the location; the additional work needed 
for points close to the boundary is small compared to that for points on the boundary or 
far away. For an integral with sources on several boundaries, the evaluation at a point 
on one surface might lead to a nearly singular case, since it could require integration 
over another surface which is close to the first. Results in [ |33| showed that the two- 
dimensional version [|3]| of the present method works well in such cases. 
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In this paper we treat single and double layer potentials for harmonic functions. The 
approach can be extended to integrals for Stokes flow as in | |28| ; see also [9. 10,2lj. This 
integration method could be applied to problems with moving interfaces. A discretiza¬ 
tion of the evolving surface must be chosen and updated. This is often done with rep¬ 
resentative marker points and triangulation. An alternative might be to use the level set 
method [j22].[26j|, representing the surface as the zero set of a function whose values are 
transported by a computed velocity. The velocity would be needed at nearby grid points, 
in order to update the level set function, rather than on the surface, and the present inte¬ 
gration method is designed to be suitable for this purpose. 

In Section[2]the quadrature rule for surface integrals is explained. In Section[3]the for¬ 
mulas for calculating single and double layer potentials are given, including the simpli¬ 
fied version for points on the surface. A brief discussion of the error estimates is included. 
Numerical examples illustrating the method with a variety of surfaces are presented in 
Section El Finally, some possible improvements for this method are discussed briefly in 
Section|5| The code that produced the examples is available on request from the first two 
authors. 


2 The quadrature method for surface integrals 

In this section we describe the computational method for integrals on implicitly defined 
closed surfaces in three space dimensions. The method is also applicable to closed curves 
in 1R 2 and to hypersurfaces in higher dimensions. Detailed proofs and extensive examples 
are given in Wilson [|30p. Our purpose is to evaluate the surface integral 

I = j T f(y)dS y (2.1) 

We assume T is a C 2 "' fl surface, m > 1, and, for now, / G C 2m (T). The method exploits 
the spectral convergence of the trapezoidal rule without the need for the user to generate 
a set of overlapping coordinate patches and associated partition of unity. Rather than 
covering T with overlapping rectangular patches, we cover T with certain overlapping 
surface sets. We define the subsets 

Tj = {xeT: |n(x)-e;| >0}, z = 1,2,3 

where n is the unit outward normal at x and {e,} 3 =1 is the standard basis for R 3 . Thus T, 
contains all points in T where n is not orthogonal to e, . 

The integration method uses a high order, patch-independent quadrature formula 
for integrals with integrands that vanish outside of a compact subset of one T,. To handle 
general integrands /, we introduce a partition of unity to find functions {/'} 3 =1 such 
that /(x) = Yji=\f l ( x ) for all x G T and /' vanishes outside a compact subset of T ( . We first 
design a universal partition of unity on the unit sphere, S = {uGR 3 :|u|= 1}. We start with 
the smooth bump function defined as b{r) = e r2 /b 2 -i) for \r \ < 1 and b(r) = 0 otherwise. 
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Next we choose a fixed angle 9 with cos 1 
we define 

w 1 {vl) = cos“ 1 (|u-e,|). 


(1/ \/3)<9<n/2. For each z'G {1,2,3} anduGS 


a l,e { u): 


b(zv l (u)/6) 

E b(ivi(u)/6) 
;=i 


Because 9 >cos“ 1 (l/ 43), the sum is always positive. Furthermore 

1. For each i = 1,2,3, we have a 1,6 G C°°(S); 

2. For all u G S, we have Ef=i C7 ' ! ' 0 ( u ) = 4 

3. For each i = 1,2,3, the function a l,e vanishes outside the compact subset 


S^ = {uGS : |u-e,| >cos0} 


To make use of the above partition of unity on the sphere for a general surface Y, we 
apply it to the unit normal n(x) G C lm (Y). The composition functions = a 1,6 on on F 
satisfy 

1. For each i = 1,2,3, we have £, l,e G C 2m (T); 

2. For all xGT, we have Ef=i4 ,0 ( x ) = 1; 

3. For each i = 1,2,3, the function C, l,e vanishes outside the compact subset of F, given 
by 

F,^ = {x G F : |n(x)-e;| >cos0} (2.2) 

Using the partition of unity, we obtain the exact formula 

3 

J T f(y ) dSy = £ J T C' e (y )/(y) dS y (2.3) 


where the integrand f'(y) = 4' 0 (y)/(y) vanishes outside the compact subset F,^ of V 
Finally, the surface integral (|2.1| can be approximated by the numerical quadrature 
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h 2 L E 

1= lx£K|i,j,s 


4' e (x)/(x) 

|n(x) -e/| 


Jz 2^ t7i ' e ( n ( x ))/( x ) 

i=lxeR h ,i,e 


n x -e,; 


(2.4) 


where 

Rh,i,e = { x £T ': |n(x)-e,| >cos0andp'( x ) £h^ 2 } 

and p l : 1 R 3 i-t 1 R 2 is the projection function defined by p 1 (x) = (X2/T3)/ p 2 (x) = (x\,X3), 
p 3 (x) = (X],X2). Rh,i,e consists of those points in F,^ that project to grid points in the 
corresponding plane. The weights l/|n(x) e,| correspond to the area elements of the 
inverse projections. 
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It is proved in [301 that I], —I = 0(h 2m ), i.e v the quadrature rule (2.4' is high order 
accurate, provided the surface Y is C 2m+1 and / is C 2m , m > 1; see Lemma 1, Theorem 2, 
pp. 9-10, and Lemma 9, Theorem 10, pp. 25-27. In effect the trapezoidal rule applies on 
coordinate patches covering each T^g. The spacing h must be small enough to resolve the 
surface. Assume T is defined as the set (p = 0 for some function cp on an open subset of 
R 3 . If Ci =min| Vcp\ and C2 = max||D 2 </>|| near the surface, we need h <h.Q = 2Cicos0/ C 2 . 
Thus if the curvature is large, h must be small. The method works if <p is known only at 
grid points; see p. 30 of |30) . 

The use of this rule requires finding the points in Rh,i,e■ Given a grid point x in a 
coordinate plane, there may be several points in T which project to x, but they are well 

pp. 


separated because of the normal condition in (2.2 1 . Consequently, as shown in 
11-18, a simple line search algorithm can be used to locate the quadrature points if h < ho- 
Briefly, to find points in R] t/3/ e, for each (j' 1 , 72 ) and 73 , check whether cp has a root x = 
{jih,j 2 h,x 3 ) with j 3 h <x 3 < (j 3 + l)h. If so, find a root. If the root is in R} h3/ g, it is unique. If 
it is not in Rh^fis reject it. In either case go to the next 73. The validity of this algorithm for 
h<hois proved in [301. 


3 Evaluation of the layer potentials 

We describe the procedure for computing a single or double layer potential at an arbitrary 
point, the most difficult case being a location off the surface but close by. For the case of 
a point on the surface we give more special versions with improved accuracy. Finally we 
discuss error estimates. 

3.1 The single layer potential 

The single layer potential on F determined by a density function ip is 

v{x) = ^G(y-x)iKy)dS y (3.1) 

where G is the fundamental solution for the Laplacian, G(x) = — 1/4tt|x|. We suppose I 
is the boundary of a bounded domain Q. To evaluate v for x close to T, we replace G with 
a smoothed, or regularized, version 

G i (y) = G(y)erf(|yj/^) = -5SML (3.2) 

where erf is the error function and G^O) = — rc^ 3 ^ 2 (25)~ 1 . The resulting error in the 
integral is 0(<5). Typically 1 < S/h < 2 . We first compute the regularized integral 

v s {x) = J^G s (y-x)ip(y)dS y . 


(3.3) 
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using the method of Sec. [2j The value obtained is not close to v(x) because of the near 
singularity. We add corrections for the regularization and discretization to improve the 
accuracy. Other regularizations could be used, but (3.2} has the advantages that it is 
simple, G$ — G decays rapidly in the far field, and manageable formulas can be found for 
the corrections described below. 

To obtain the corrections we first find z G T, the closest point on the surface to x, and 
set x = z + bn. Here n is the unit outward normal to T at z; b < 0 if x G Q and b > 0 if x G Q c 
(the complement of the closure of Q in R 3 ). The correction for regularization of the single 
layer potential is 


71 = - (l+HAS) ip(z) 


| A | erfc | A | — 


„-A 2 


V^J 


(3.4) 


Here erfc(r)=l—erf(r), A =b/5 and H is the mean curvature at z, H=(k\+K 2 ) /2, where K\ 
and K 2 are the principal curvatures. Formulas for computing needed geometric quantities 
such as H are given in Appendix B, and the sign convention for H is explained. 

The discretization correction is a rapidly convergent infinite sum resulting from the 
Poisson summation formula, applied in each T„ i = 1,2,3. Let 


E(p,q)=e 2,,q erfc(p+q)+e 2pq erfc( — p+q). (3.5) 


Also let Q = {n = ( n\,ti2 ) G Z 2 : H2 > 0 or (ti2 = 0 andni > 0 )}. Now suppose z lies in a 
system of coordinates, say oc = (ai,«2); in our case, zGTy for one or more of k= 1,2,3, and 
the coordinates at z are p k ( z) G R 2 . Let g l i be the inverse metric tensor at z, and for n GQ, 
define \\n \\ 2 = Y t i / jg iq nitij. Also write the coordinates of z as (^1,0:2) = (mi r m.2)h+(v\ r V2)h 
where mi,m2 are integers and 0 < V\,V2 < 1 . Here (1/1A2) = v ' k) arid ||n|| depend on the 
choice of k. The discretization correction is 


72 = ^( z )E J^C k '%z)cos{2nn-vW)^E(A,n5\\n\\/h) 


(3.6) 


k=lneQ 


with [, k ' e as in Sec. [ 2 ] Finally, the computed value v of the integral (3.11 is v ~ Vs -V ’/~i + 'Ti, 
where v$ is the value of (3.3} obtained by the quadrature rule. 


3.2 The double layer potential 

The double layer potential has the form 

^)=j r dG( i < p( y ) ds y ( 3 - 7 ) 

It is discontinuous at F. If x is close to F, we find the closest point z and distance b as 
before. We use Green's identities to reduce the singularity and then regularize the kernel, 
obtaining 


^ (X) = J T dGS ^W(y)-?( z )]d5y+X^) 


( 3 . 8 ) 
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Here x = 1 f° r x £ H, x = 0 on Q c , x = \ on T. To form / 3 n we use 

smooth function Gs introduced in \3.2\, 

the gradient of the 

VG,(y) = VG( y ) S (| y |/,i) = i ^3 S (y|/<5) 

(3.9) 

with 


s(r)=erf(r) - -=re r ~ 

\ n 

(3.10) 

We compute the integral in (]3.8|> as in Sec. |2land again add corrections. The regularization 

correction for (|3.8|> is 


M-<5 2 (A s <p) 4 ||A|erfc|A| . 

(3.11) 


Here Ag cp is the surface Laplacian of cp at z, which is expressed in coordinates as 



with g = det gij. The discretization correction is similar to 7? but involves dcp/da r , r = 1,2, 
the coordinate derivatives of cp evaluated at z. In our case, on T^g, are the two 

components of x= [xi,X 2 ,x^) other than x k , and we write these derivatives as d^cp(z). 
The correction is 

(3.i2) 

Z k=lr=l 


where 


4^ = ^sin(27rn-v^)^j-^pE(A,7r£||n||//i) 

ngQs=l W n W 


with ||n|| as before. The computed value of (3.7' is w^uvs+Mi+Mi- 


(3.13) 


3.3 The potentials evaluated on the surface 

We now treat the important special case of evaluation at a point x on the surface T. For 
this case, in contrast to the nearly singular case, it is not difficult to modify the regularized 
kernels to have higher accuracy by imposing moment conditions. Thus no corrections are 
needed for regularization. The method is easier to use than in the general case, and the 
error is typically smaller, as seen in the examples in Section 4. A strategy for producing 
these improved kernels from the ones already chosen in ( |3.2| and ( |3.9| is described in (2| 
for the single layer and in Q for the double layer. For G$ as defined in |3.2| , the error in 
regularization, i.e. the difference between the integrals in ( |3.3| and ( |3.1| , is O(S). For the 
new version of Gs the error is 0(S 5 ) for the special case of evaluation on the surface, and 
similarly for the double layer. 
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To evaluate the single layer potential 0 at x G T, we use the new version of G$ with 
0 (^ 5 ) accuracy. 


G s {y) = - S( j yl /f ) / s(r)=erf(r) + -^(5r-2r 3 )e - r2 


4/r | y | 


3 sfn 


(3.14) 


and G tf (0) = — (4/3)71 3/2 £ 1 . In place of the corrections (3.4 3.6 1 we have 7i =0 and 

3 


T 2 = ^-ip(z) I] ^(z)cos( 2 /rn-v (,c) )F(^) 
71 k=lneQ 


where, with £ = 27z\\n\\6/h and ||n|| as before. 


F(£) = ^erfc(£/2) + 7 r 1/2 ^e ^ /4 (^l + ^- 


l 


(3.15) 


(3.16) 


The derivation of (3.14' is similar to that in Sec. 2 of ||2) for an 0(5 3 ) version, with the 
formula ( 3.15} corresponding to (3.28), (2.2 3) in |4j. 

To evaluate the double layer potential (3.7 1 at x G T, we use (3.8 1 with X = 2 ant l VC d 
of the form ( ]3.9| ) but with ( |3 .10 1 replaced by 


s(r) =erf(r) —( r-%-\e 




3 J 


(3.17) 


In this case no corrections are needed, that is, = J\f 2 = 0. Formula ( 3.17} was derived 
in |4||, p. 607. 


3.4 Error analysis 

For the general case of points close to T, the corrections M\, N 2 for the double layer 
were derived in |4j. Those for the single layer can be found similarly; we include a brief 
derivation of T\ in Appendix A. The discretization corrections are based on the Poisson 
Summation Formula. After applying both corrections to either the single or double layer 
potential, the remaining error has the form (cf. Theorem 1.2 in |4j) 

£<C 1 5^ + C 2 h 2 e- Co(s,h)2 (3.18) 

as 5, h —> 0, assuming 5/h is bounded below, with Ci,C 2 depending on derivatives of 
the surface and density functions. The two terms represent the regularization error and 
discretization error, respectively. The constant Co is determined by the choice of local 
coordinate system. It is important for accuracy that Co does not become small, so that the 
sums in the discretization corrections converge rapidly, and so that the second error term 


c 0 < 7 = minW'Tih 

l fc l =1 i,j 


in (3.181 is comparable to the first in practice. It was shown in [4J, Sec. 3, that the estimate 
(3.181 holds provided 
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(A factor 1/2 in Cq 
With coordinates (ai 




is the inverse metric tensor. 


in [41 was arbitrary.) Here g'i = (gij)~ 

i],tt2), gij = T, ■ Tj, where T) = dx/dttj, j = 1,2 are the tangent vectors. 
If the coordinate system distorts distances significantly, 7 could be small, and thus the 
accuracy of the method depends on the choice of coordinates. 

In our case, the coordinate systems are those determined by the projections. For [" 3 , 
the coordinates are Xi,* 2 , with x^=f{x\,X 2 ), and similarly for [/, Y 2 . From the expression 
for g 1 -' in Appendix B, it is not difficult to see that in this cas e 7 = (1 +/|) ~ ' /2 — |n ■ ea | > 


cost? in Yy), where 9 is the angle chosen in Sec. |2| Thus (3.181 holds with Co ~ n I cos I 6, 
and the exponential in (3.181 can be made quite small. For example, if 5 = 2 h and 9 = 60°, 
the exponential is .00005; if S=2.5h and 0=70° it is .0007. In practice the accuracy is about 
0(/z 3 ) for usual values of h, with proper choice of parameters. Alternatively, we could 
take 5 = CM for any q < 1 and thereby obtain convergence as h —>0 with order 0 (/j 3<? ). 

For the sums in (3.6 1 , (3.131 we only need a few terms because of the rapid decay as 
n increases. In the corrections we may evaluate H and A$cp not at the closest point z 
but rather at a neighboring grid point, using formulas in Appendix B. The 0(h ) errors in 
these quantities do not change the order of accuracy of the corrections. Similarly, in the 
discretization corrections, the g*i and d? <P only need to be computed within 0(h). 

For the special case with x G Y, the first term in the error estimate (3.181 improves to 
S 5 . (See Theorem 1.1 of Q for the double layer.) Thus we can take S/h somewhat larger 
to reduce the discretization error in the second term. If 5 = CM, q < 1, then the error is 
0(h 5t i) as q— > 0 . 


4 Numerical Results 

This section presents examples evaluating the sum of a double layer and single layer 
potential on five different surfaces, 

u(x) = w(x)+v{x) = J^ dG ^ X \ (y)ds y -J^G(y,x)xp(y)d.Sy. 

In all the examples the potential u is chosen to be 

, . I (sinx+sim/)e z if (jc,y,z) G O 

u(x,y,z) = < 

[O if (x,y,z) GO c 

The densities cp and ip are determined by the jumps in u and du/dn. The integrals are 
calculated given these densities, and the result is compared with the exact u. This choice 
of test problem allows us to have an exact solution with an arbitrary surface. 

In each example, the domain Q is embedded into a cubic box B=(—L,L) 3 with L=l.l. 
The box is partitioned into a uniform grid Tj t with mesh parameter h = 2 L / N, the width 
of a grid cell. We call a grid node irregular if the stencil of the second-order Laplacian A ; 2 
crosses the boundary F; otherwise it is regular. 
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The numerical values of u(x) = zv(x) + i’(xj are first computed at the irregular grid 
nodes, as well as the neighboring nodes in their stencil, using the procedure of Sec. [3] 
The sums Wg, Vg for the smoothed potentials are found and the corrections are added. 
The summation is done using a slight modification of the treecode algorithm of Duan 
and Krasny 0 , which was designed for use with Ewald summation. (The kernel in 
their code has a factor of erfc rather than erf, and their solutions are periodic rather than 
in free space.) In the treecode we chose the degree of Taylor polynomials p = 12, the 
separation parameter s = .5, and the capacity, or maximum number of points in a leaf, 
N 0 = 20. 

Having calculated w(x) at grid nodes x near T, we can now find values at all the 
regular grid nodes of 7 1, by inverting the discrete Laplacian, using a procedure suggested 
in |20| . Let u^g denote the value of u already computed at the nodes close to F as nearly 
singular integrals. We formulate a Poisson problem for an approximation Uj, on T to the 
exact u, 

J A;, U] h g at irregular grid nodes 

= < ' 

10 at regular grid nodes 

We solve for uj, with a fast Poisson solver on the box B with zero boundary condition. 
At the regular nodes, the truncation error is 0(h 2 ), since the exact solution is a smooth, 
harmonic function away from F. For the irregular nodes, Uj h g — u is about 0(/; 3 ), so that 
^h u h,s~ A;,m = 0(h 5 /h 2 ) = 0(h). In summary, A/ ; (m/,— u) is 0(/z 2 ) at regular nodes and 
0(h) at irregular nodes. Despite the first order truncation error near F, the resulting «/; is 
second order accurate, i.e., uj l —u=0(h 2 ), uniformly on T, a fact proved in |5]. While this 
method is efficient, more accurate values could be computed as integrals or otherwise. 

We also computed the values of u (x) at grid nodes x on F using the method in Section 
3.3. (The exact value is the average of the inside and outside limits.) In this case we used 
direct summation to provide an unambiguous test of the accuracy. 

For Example 1, the surface T is an ellipsoid given by 


x 2 y 2 z 2 

+ 7T + =T = 1 




b 2 


with a = 1, b = . 8 , c = . 6 , and rotated by an orthogonal matrix to test the effect of grid align¬ 
ment. Example 2 is a thinner ellipsoid, with a = 1, b = c = .4, without rotation. Example 3 
is a torus 


^x 2 +y 2 - C ) 2 +z 2 = a 2 

with a = .3 and c = .7. Example 4 is a molecular surface with four atoms, similar to one of 
the definitions in ( 8 ), 

£]exp(|x —x fc | 2 /r 2 ) = c 
k =1 
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Figure 1: The rotated (1,.8..6) ellipsoid 


Here the centers are (y/3/3,0,— x/6/12), (—\/3/6,.5,—\/6/12), (—\/3/6,—.5, —\/6/12), 
(0,0,\/6/4) and r = .5, c = .6. Example 5 is a surface obtained by revolving a Cassini oval, 

(x 2 +y 2 +z 2 +a 2 ) 2 - 4a 2 (x 2 +y 2 ) = b 4 


with a = .65 and b = .7. 

The errors for each example are presented in the tables with N = 64, 128 or 256 and 
with S/h = 1, 2, or 3. Both L 2 and maximum, or L°°, errors are given. They are displayed 
first for the irregular grid points, then for the regular grid points, and finally for the 
quadrature nodes on the surface. In the L 2 norms, points are given equal weight, and 
thus for the quadrature nodes this measure effectively gives extra weight to the overlap 
regions. In all the examples the angle 6 in the partition of unity on the sphere is 70°. The 
final table gives the number of quadrature nodes on each surface with N = 256. 

For the irregular points the smallest errors are generally for S/h = 1, while accuracy 
approaching 0(h 3 ) is observable with S/h = 2 or 3. (We have found that errors are larger 
with S/h < 1.) The errors at the regular points is 0(// 2 ) as expected. The errors at the 
quadrature points are generally smaller than those at the irregular points for S/h = 2 or 
3, but not for S/h = 1; this reflects the fact that the method of Section 3.3 for evaluation on 
the surface improves the smoothing error directly but the discretization error is improved 
indirectly by the smoothing. With S/h = 3 the errors on the surface decrease rapidly 
with refinement. We repeated the computations on the surface with angle 60° and found 
similar but slightly larger errors. To test further refinement, we computed u on the surface 
for Example 1 with N=512. With angle 70° and S/h= 2, the L 2 and L 00 errors were 3.03E-8 
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Figure 3: The torus 
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Figure 5: The Cassini oval surface 
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and 7.43E-7; with S/h =3 they were 3.09E-10 and 1.33E-8. With angle 60", they were 6.83E- 
10, 3.82E-8 for 5/h = 2 and 2.42E-10, 1.09E-8 for 5/h = 3. We conclude that for practical 
use we can reliably choose S/h = 2 for values at points near the surface and 6 / h = 3 with 
the special method of Sec. 3.3 for points on the surface. 

In these calculations we used analytical values of g t] since the surface was specified. 
However these could easily be computed from grid values of a level set function (see 
Appendix B) and are needed only to 0(h) accuracy. We used analytical values of the 
densities cp and ip at the quadrature points and computed other values numerically from 
these. We did so because in solving an integral equation we would only know values 

rjA 

at the quadrature points. Occasionally the computation of d) cp in (3.12i failed for lack 
of nearby quadrature points near the edge of the support of one In such a case we 
set the contribution to zero since £ k,e is very small there. We also treated the ellipsoid of 
Example 1 without rotation; the errors were similar to those displayed with the rotation 
but slightly smaller. We tried an ellipsoid thinner than in Example 2; we found that the 
interpolation stencil needed for the corrections T\, N\ failed with N =64 but worked with 
larger N. 


Table 1: Errors for the rotated (1,.8,.6) ellipsoid. In each table the angle 8 is 70° and h = 2.2/N for the N 3 
grid. L 2 and L°° errors are displayed for (i) irregular grid points near the surface, (ii) regular grid points, and 
(iii) quadrature points on the surface. 


5 

grid 

lle irres lh 
W e h 112 

lh meg l 

W C h 11“ 

lle reg IU 

II 112 

IIHr 11“ 

11 c l ua d 11 

W e h 112 

11 c l uac ^ 11 

II e h 11“ 

h 

64 3 

3.32E-5 

2.57E-4 

3.29E-5 

5.05E-4 

9.80E-5 

1.24E-3 

128 3 

4.14E-6 

3.54E-5 

7.96E-6 

1.07E-4 

2.10E-5 

3.33E-4 

256 3 

9.91E-7 

6.55E-6 

2.10E-6 

2.92E-5 

5.36E-6 

8.69E-5 

2 h 

64 3 

1.39E-4 

8.64E-4 

1.82E-4 

2.39E-3 

3.15E-5 

4.32E-4 

128 J 

1.78E-5 

1.14E-4 

4.35E-5 

4.67E-4 

1.71E-6 

3.94E-5 

256 3 

3.33E-6 

1.39E-5 

1.13E-5 

1.36E-4 

1.38E-7 

4.69E-6 

3 h 

64 3 

4.78E-4 

2.91E-3 

3.65E-4 

5.06E-3 

2.59E-5 

2.83E-4 

128 3 

5.75E-5 

3.82E-4 

8.79E-5 

9.45E-4 

1.13E-6 

1.68E-5 

256 3 

7.84E-6 

4.78E-5 

2.27E-5 

2.70E-4 

1.67E-8 

6.36E-7 


5 Discussion 

The numerical results illustrate the performance of the method and are in general agree¬ 
ment with the qualitative predictions. They show that reasonable accuracy can be ob¬ 
tained with moderate resolution, and the observed order of accuracy gives confidence 
that the errors will reduce with further refinement. Here we comment on possible im¬ 
provements. 

As noted in Section 3, the discretization error depends on the angle 9 in the partition 
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Table 2: Errors for the (1,.4,.4) ellipsoid 


5 

grid 

lle irres lh 
In/; l|2 

lh mes ll 

II e h II 00 

K s h 

IKloo 

11 c l uai ^ 11 

II e h Il2 

11 °l ua( ^ 11 

II 6 /, ll°° 

h 

64 3 

4.46E-5 

3.27E-4 

2.62E-5 

4.39E-4 

1.48E-4 

9.26E-4 

128 3 

8.53E-6 

8.73E-5 

6.79E-6 

1.21E-4 

3.25E-5 

3.07E-4 

256 3 

1.08E-6 

1.39E-5 

1.69E-6 

3.15E-5 

7.35E-6 

6.70E-5 

2 h 

64 3 

3.29E-4 

2.33E-3 

1.69E-4 

2.05E-3 

4.95E-5 

2.66E-4 

128 3 

4.29E-5 

2.73E-4 

4.07E-5 

4.17E-4 

6.82E-6 

9.36E-5 

256* 

5.11E-6 

3.57E-5 

1.04E-5 

1.05E-4 

4.04E-7 

1.00E-5 

3 h 

64* 

1.12E-3 

6.85E-3 

3.37E-4 

4.29E-3 

4.59E-5 

2.35E-4 

128 3 

1.43E-4 

9.62E-4 

8.07E-5 

7.85E-4 

5.22E-6 

5.20E-5 

256 3 

1.76E-5 

1.23E-4 

2.09E-5 

2.22E-4 

2.03E-7 

3.91E-6 


Table 3: Errors for the torus 


<5 

grid 

lle irres lh 

II e h 112 

11 p me S 11 

ll e /i H°° 

|| e res |L 
In/, 11 2 

iknioo 

11 4 uac * 11 

ll e /, l|2 

lid- 

h 

64 3 

7.19E-5 

3.57E-4 

5.17E-5 

5.03E-4 

1.48E-4 

1.29E-3 

128 3 

8.61E-6 

7.56E-5 

9.04E-6 

1.14E-4 

3.14E-5 

3.05E-4 

256* 

1.02E-6 

9.61E-6 

2.16E-6 

2.01E-5 

7.02E-6 

7.50E-5 

2 h 

64* 

2.42E-4 

7.94E-4 

1.74E-4 

1.42E-3 

8.08E-5 

4.16E-4 

128* 

2.89E-5 

9.54E-5 

4.41E-5 

3.31E-4 

8.53E-6 

8.34E-5 

256* 

3.52E-6 

1.25E-5 

1.12E-5 

8.85E-5 

4.76E-7 

7.46E-6 

3 h 

64* 

8.17E-4 

2.68E-3 

3.17E-4 

3.01E-3 

6.35E-5 

2.80E-4 

128* 

9.92E-5 

3.28E-4 

8.51E-5 

6.93E-4 

7.05E-6 

4.85E-5 

256* 

1.22E-5 

4.22E-5 

2.18E-5 

1.84E-4 

2.46E-7 

2.80E-6 


of unity on the unit sphere, defined in Sec. 2. We need 9 > 55° to cover the sphere. As 
6 increases toward 90° we expect the accuracy to deteriorate because of the dependence 
of the discretization error on 6, as explained in Sec. 3.4. Here we used 6 = 70° as a 
compromise between the extremes. In our experiments the errors were not very sensitive 
to the choice of angle for 60° <6 < 80°. We found slightly larger errors with 9 = 60° than 
for 70°. A possible explanation is that the gradient of the partition of unity functions is 
larger for the smaller angle. It is unclear whether this can be improved or is an inherent 
limitation. 

In the discretization corrections T 2, J\f2 we summed over n = [n\, H 2 ) with \tij\ < 
20, but the number of terms actually needed is much smaller. In fact for 5/h> 2 these 
corrections are usually negligible. They could be modified to include an estimate of the 
number of terms needed to avoid unnecessary work. 

The treecode of flTfl cannot be used directly with the method of Section 3.3 for eval¬ 
uation on the surface because of the differences in the regularized kernels. However, 
the treecode could be applied in this case by modifying the recurrence formulas for the 
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Table 4: Errors for the molecular surface 
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grid 

||g irres |L 

W c h 112 

iiciioo 

||g res |L 

W e h 112 

lh res ll 

W e h 11“ 

11 4 uac * 11 

\\ e h Il2 

ikriu 
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64 3 

6.84E-5 

4.16E-4 

4.02E-5 

4.81E-4 

1.61E-4 

1.47E-3 

128 3 

5.98E-6 

5.55E-5 

1.09E-5 

1.50E-4 

3.24E-5 

3.24E-4 

256 3 

1.03E-6 

1.30E-5 

2.73E-6 

3.96E-5 

7.60E-6 

8.33E-5 

2 h 

64 3 

3.30E-4 

1.53E-3 

2.51E-4 

3.73E-3 

6.80E-5 

6.07E-4 

128 3 

3.99E-5 

1.81E-4 

6.02E-5 

7.37E-4 

3.43E-6 

6.78E-5 

256 3 

4.96E-6 

2.36E-5 

1.45E-5 

1.57E-4 

2.32E-7 

4.34E-6 

3 h 

64 3 

1.11E-3 

5.12E-3 

5.08E-4 

8.56E-3 

6.35E-5 

4.35E-4 

128 3 

1.38E-4 

6.26E-4 

1.18E-4 

1.58E-3 

2.01E-6 

3.46E-5 

256 3 

1.72E-5 

8.17E-5 

2.80E-5 

3.27E-4 

5.40E-8 

1.24E-6 


Table 5: Errors for the Cassini oval surface 
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grid 
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1 1 p 1TIe S 1 1 

II e h 11“ 

\Kh 
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11 ( l ua< ^ 11 
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11 c l ua d 11 

II e h 11“ 
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64 3 

4.87E-5 

2.94E-4 

3.40E-5 

3.60E-4 

1.20E-4 

1.01E-3 

128 3 

3.78E-6 

3.07E-5 

7.25E-6 

6.59E-5 

2.37E-5 

2.27E-4 

256 3 

6.82E-7 

5.75E-6 

1.82E-6 

1.68E-5 

5.62E-6 

5.75E-5 

2 h 

64 3 

2.02E-4 

8.64E-4 

1.55E-4 

1.20E-3 

5.18E-5 

3.53E-4 

128 3 

2.46E-5 

1.20E-4 

3.80E-5 

3.38E-4 

3.15E-6 

3.76E-5 

256 3 

3.10E-6 

1.56E-5 

9.59E-6 

8.17E-5 

2.19E-7 

3.63E-6 

3 h 

64 3 

6.83E-4 

2.65E-3 

2.97E-4 

2.58E-3 

4.47E-5 

2.20E-4 

128 3 

8.61E-5 

3.91E-4 

7.35E-5 

7.50E-4 

1.84E-6 

1.86E-5 

256 3 

1.08E-5 

5.13E-5 

1.86E-5 

1.78E-4 

4.65E-8 

6.96E-7 


Taylor coefficients as derived in [111. This could be done in future work. 


Other than [ 11 ], fast summation methods have not been developed specifically for 


regularized kernels. Among existing codes, one possible alternative that might be used 
in the present computations is the kernel-independent fast multipole method (KIFMM) 


of L. Ying, G. Biros, D. Zorin [31]. We have calculated examples using this code, even 
though it was not intended for regularized kernels. We found difficulty maintaining 
good accuracy, especially with larger 5/h, perhaps because the regularization degrades 
the accuracy of the linear problems solved in the KIFMM. We emphasize that this is a 
use of the KIFMM for which it was not intended. A summation method of fast multipole 
type designed particularly for these regularized kernels could improve the efficiency of 
this method without loss of accuracy. 
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Table 6: The number of quadrature nodes with N = 256 


Example 

1 

2 

3 

4 

5 

Number of nodes 

144388 

70790 

142168 

126789 

133014 
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A Regularization error for the single layer potential 


The regularization correction for the single layer potential, evaluated at a point near the 
surface, can be derived using the method of |4], Sec. 2. With G$ and G as in ( |3.2| , we 
approximate the error in the single layer potential with density xp, evaluated at a point 
x near the surface. Since the error is local, we write it as an integral in one coordinate 
patch, regarding xp as a function of a = (^ 1 ,^ 2 ), 

e = J [G s (y(a)-x)-G(y(a)-\)}xp(a)dS(a). (A.l) 

For simplicity, we will assume that x is along the normal line from 0 G T, so that x = brio 
for some b, where no is the unit normal at 0. We also assume the coordinates are chosen 
so that a(0) =0, gij{oc ) = <G+0(|a| 2 ), and the tangent vectors T\, T2 have the directions of 
principal curvature. Thus 


e = 


' erfc (r/d) 


47T, 


xp(a)dS(a), r=\y(a) — x\ 


(A.2) 


Proceeding as before, we make a near-identity coordinate change a —> £ such that 
|£| 2 +b 2 = r 2 . We get 


€= ■ 


47T . 


' erfc(y / |£| 2 + fr 2 / 5) 

( l £| 2 + & 2 ) 1/2 


zv{£,b)dcp, 


with 


w(g,b) = ip 


da 

dC 


I Ti x T 2 1 


(A.3) 


(A.4) 


We will see that we can neglect terms in w of the form 0(|£| 2 + fc> 2 ). We can approximate 
xp in £, with leading term xpo = xp(0), 


xp = xp 0 +xp j (l+bcj/2^ j + 1 1 xp ij ^j+Om+O(b 3 ). (A.5) 
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where (7 = 7 Ci £ 2 /|£| 2 + k : 2 £ 2 /|£| 2 and Ki,jc 2 are the principal curvatures at 0. For the other 
two factors in w, we have 


det(9a/9£) = l + b( 7 + 0 (|£| 2 )+ 0 (b 2 ). (A. 6 ) 


and 

|TixT 2 |=1+0(|^| 2 ). (A. 7) 

In the ^-integral for e, the odd part of iv will contribute zero. Thus we can replace w 
with an approximation to its even part. Combining the three factors above, we get 

w™ n (Z,b) = xp 0 (l + bci)+O(\Z\ 2 + b 2 ) (A. 8 ) 


We now substitute in the integral, change to polar coordinates, and substitute |£| =5s and 
b = 5 A to obtain 

e= S -Ul +«H) £ er *0^^ sis+O(?) (A.9) 

with H = (k\ +k 2 ) /2, the mean curvature. With r = \/s 2 + A 2 , and sds = rdr, the integral 
simplifies to 


r Crisis 

Jo r 



= e ^ / \frc — | A | erfc | A | 


and finally 

e = ^ 0(1 +SAH) (e~ x2 / y/n- |A|erfc|A|) +0(£ 3 ) 
leading to the correction (|3.4|>. 


(A.10) 


(A.11) 


B Formulas for Monge Patches 


We summarize formulas needed for the corrections of Sec. [3]when applied in a coordinate 
system such as *3 =f(x\,x 2 ), often called a Monge patch. Given /, let /, = df/()x u i = 1,2, 
and similarly for a second derivative fq. The metric tensor (gq) and its inverse (g 1 - 1 ) are 


(gij) = 


1+/1 2 / 1/2 A 

/ 1/2 1 +fl)' 


(s li ) = 


1 

g 


i+/l 

-/ 1/2 


-/ 1/2 \ 

1+/1 2 y 


where g = det (gq) = 1 +/ 2 +/ 2 • The Gauss curvature is 


K=g 2 (/ii/22-/i 2 2 ) 


The mean curvature is 


H ~±lg 3/2 [( 1 +/ 2 )/ll + (l+/l)/ 22 - 2 / 1 / 2 / 12 ] 
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where the sign is + if x$> f(xi,X 2 ) outside and — otherwise. 
The surface Laplacian has the general formula 




A s<P = E^= 3 ; =Ls 11 (?&&) +E c ^' < P 




where 



With some calculation we find 




and subsequently 


Cl =gAi [2/1/2/12 - (1+/D/11 - (1 +fi)fn] = T— Hf, 


Suppose the surface is defined by (p(xi,X2,X3) =0, with (p >0 outside. Near a given 
point, there is at least one Monge patch; suppose we can solve for X3 =f(x\,x 2) as above. 
By differentiating implicitly we get /, = —</>, / t/y etc. We can use these to express gjj and 
g tJ . A more convenient expression for the mean curvature H on a surface {0=0} is based 
on the classical formula 

2H = -V-n = -V-(V0/|V0|) 

If we carry out the differentiation we get 

2H = -\V(p\~ 3 (cpatf -cpicpjcpi^ 

summed over i,j. After canceling and combining terms, we obtain a formula such as 
in(2g,p. 12. 


References 

[1] K. E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cam¬ 
bridge University Press, 1997. 

[2] J. T. Beale, A convergent boundary integral method for three-dimensional water waves. 
Math. Comp., 70 (2001), 977-1029. 

[3] J. T. Beale and M. C. Lai, A method for computing nearly singular integrals, SIAM J. Numer. 
Anal., 38 (2001), 1902-1925. 

[4] J. T. Beale, A grid-based boundary integral method for elliptic problems in three-dimensions, 
SIAM J. Numer. Anal., 42 (2004), 599-620. 

[5] J. T. Beale and A. T. Layton, On the accuracy of finite difference methods for elliptic problems 
with interfaces, Commun. Appl. Math. Comput. Sci., 1 (2006), 91-119. 


21 


[6] J. Bremer and Z. Gimbutas, A Nystrom method for weakly singular integral operators on 
surfaces, J. Comput. Phys., 231 (2012), 4885-4903. 

[7] O. Bruno and L. Kunyansky, A fast, high-order algorithm for the solution of surface scatter¬ 
ing problems: Basic implementation, tests, and applications, J. Comput. Phys., 169 (2001), 
80-110. 

[8] M. Chen and B. Lu, TMSmesh: A robust method for molecular surface mesh generation 
using a trace technique, J. Chem. Theory Comput., 7 (2011), 203-12. 

[9] R. Cortez, The method of regularized Stokeslets, SIAM J. Sci. Comput., 23 (2001), 1204-25. 

[10] R. Cortez, L. Fauci and A. Medovikov, The method of regularized Stokeslets in three dimen¬ 
sions: analysis, validation, and application to helical swimming, Phys. Fluids, 17 (2005), 
1-14. 

[11] Z.-H. Duan and R. Krasny, An Ewald summation based multipole method, J. Chem. Phys. 
113, (2000), 3492-5. 

[12] M. Ganesh and I. G. Graham, A high-order algorithm for obstacle scattering in three dimen¬ 
sions, J. Comput. Phys., 198 (2004), 211-42. 

[13] I. G. Graham and I. H. Sloan, Fully discrete spectral boundary integral methods for 
Helmholtz problems on smooth closed surfaces in IR 5 , Numerische Mathematik, 92 (2002), 
289-323. 

[14] W. Hackbusch, Integral Equations, Theory and Numerical Treatment, Birkhauser, 1995. 

[15] J. Helsing and R. Ojala, On the evaluation of layer potentials close to their sources, J. Com¬ 
put. Phys., 227 (2008), 2899-2921. 

[16] J. Helsing, A higher-order singularity subtraction technique for the discretization of singular 
integral operators on curved surfaces, preprint, 2013. 

[17] A. Klockner, A. Barnett, L. Greengard and M. O'Neil, Quadrature by expansion: A new 
method for the evaluation of layer potentials, J. Comput. Phys., 252 (2013), 332-349. 

[18] K. Lindsay and R. Krasny, A particle method and adaptive treecode for vortex sheet motion 
in three-dimensional flow, J. Comput. Phys., 172 (2001), 879-907. 

[19] O. Marin, O. Runborg and A.-K. Tornberg, Corrected trapezoidal rules for a class of singular 
functions, IMA J. Numer. Anal., 34 (2014), 1509-1540. 

[20] A. Mayo, Fast high order accurate solution of Laplace's equation on irregular regions, SIAM 
J. Sc. Statist. Comput., 6 (1985), 144-157. 

[21] H.-N. Nguyen and R. Cortez, Reduction of the regularization error of the method of regu¬ 
larized Stokeslets for a rigid object immersed in a three-dimensional Stokes flow, Commun. 
Comput. Phys., 15 (2014), 126-152. 

[22] S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2003. 

[23] C. Pozrikidis, Interfacial dynamics for Stokes flow, J. Comput. Phys., 169 (2001), 250-301. 

[24] C. Pozrikidis, A Practical Guide to Boundary Element Methods with the Software Library 
BEMLIB, C.R.C., 2002. 

[25] S. Sauter and C. Schwab, Boundary Element Methods, Springer, 2010. 

[26] J. Sethian, Level Set Methods and Fast Marching Methods, Cambridge Univ. Press, 1998. 

[27] M. Taus, G. Rodin and T. J. R. Hughes Isogeometric analysis of boundary integral equations. 
Math. Models Methods Appl. Sci. 26 (2016), 1447-80. 

[28] S. Tlupova and J. T. Beale, Nearly singular integrals in 3d Stokes flow, Commun. Comput. 
Phys., 14 (2013), 1207-27. 

[29] S. K. Veerapaneni, A. Rahimian, G. Biros and D. Zorin, A fast algorithm for simulating vesi¬ 
cle flows in three dimensions, J. Comput. Phys., 230 (2011), 5610-34. 

[30] J. R. Wilson, On computing smooth, singular and nearly singular integrals on implicitly 



22 


defined surfaces, Ph.D. thesis, Duke University (2010), 
http://search.proquest.com/docview/744476497 

[31] L. Ying, G. Biros and D. Zorin, A kernel-independent adaptive fast multipole algorithm in 
two and three dimensions, J. Comput. Phys., 196 (2004), 591-626. 

[32] L. Ying, G. Biros and D. Zorin, A high-order 3d boundary integral equation solver for elliptic 
pdes in smooth domains, J. Comput. Phys., 219 (2006), 247-75. 

[33] W.-J. Ying and J. T. Beale, A fast accurate boundary integral method for potentials on closely 
packed cells, Commun. Comput. Phys., 14 (2013), 1073-93. 

[34] W.-J. Ying and W.-C. Wang, A kernel-free boundary integral method for implicitly defined 
surfaces, J. Comput. Phys., 252 (2013), 606-624. 

[35] W.-J. Ying and W.-C. Wang, A kernel-free boundary integral method for variable coefficients 
elliptic PDEs, Commun. Comput. Phys., 15 (2014), 1108-1140. 

[36] A. Z. Zinchenko, M. A. Rother and R. H. Davis, Cusping, capture, and breakup of interacting 
drops by a curvatureless boundary-integral algorithm, J. Fluid Mech., 391 (1999), 249-92. 



